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Abstract. The goal of this paper is to design optimal multilevel solvers for the finite 
element approximation of second order linear elliptic problems with piecewise constant 
coefficients on bisection grids. Local multigrid and BPX preconditioners are constructed 
based on local smoothing only at the newest vertices and their immediate neighbors. 
The analysis of eigenvalue distributions for these local multilevel preconditioned sys- 
tems shows that there are only a fixed number of eigenvalues which are deteriorated by 
the large jump. The remaining eigenvalues are bounded uniformly with respect to the 
coefficients and the meshsize. Therefore, the resulting preconditioned conjugate gra- 
dient algorithm will converge with an asymptotic rate independent of the coefficients 
and logarithmically with respect to the meshsize. As a result, the overall computational 
complexity is nearly optimal. 



1. Introduction 

In this article, we construct robust multilevel preconditioners for the finite element 
discretization of second order linear elliptic equations with strongly discontinuous coef- 
ficients. We extend corresponding results on uniform grids [|55l to locally refined grids 
obtained by bisection methods. We consider the following model problem : 

r -V-(aVw) = /inl], 

\ u = goon Yd, = on 

where f2 G M'^ is a polygon (for d = 2) or polyhedron (for d = 3) with Dirichlet 
boundary Yd and Neumann boundary Fjv such that U Fat = dVl. The diffusion 
coefficient a = a{x) is piecewise constant. More precisely, the domain is partitioned 
into M open disjoint polygonal or polyhedral regions f2j (i = 1, ■ ■ ■ , M) and 

a\n, = «i, i = 1,. . . ,M 

where each is a positive constant. The regions {i = 1, ■ ■ ■ M) may possibly have 
complicated geometry but we assume that they are completely resolved by an initial 
triangulation %■ Our analysis can be carried through to more general cases when a{x) 
varies moderately in each subdomain and to other types of boundary conditions in a 
straightforward way. 

The problem ( |1.1[ ) belongs to the class of interface problems or transmission problems, 
which are relevant to many applications such as groundwater flow [|29l , electromagnet- 
ics [|27l . semiconductor modeling ll22ll3T1l . and fuelcells ll48l . The coefficients in these 
applications may have large jumps across interfaces between regions with different ma- 
terial properties, i.e. J(a) := maxj aj/ minj aj ^ 1. Due to J(a) and the mesh size, the 
finite element discretization of ( |1.1[ ) is usually very ill-conditioned, which leads to de- 
terioration in the rate of convergence of multilevel and domain decomposition methods 
[I3l|26l|45l. 
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In some special situations, one is able to show the (nearly) uniform convergence of the 
multilevel and (overlapping) domain decomposition methods (see [fT2ll46ll47ll23l[36l for 
examples). For general cases, one usually need some special techniques to obtain robust 
iterative methods, (cf. [16lllQl|251[r]). Recently in [55l|59l|, we analyzed the eigenvalue 
distributions of the standard multilevel and overlapping domain decomposition precondi- 
tioned systems, and showed that there are only a small fixed number of eigenvalues that 
may deteriorate due to the discontinuous jump or mesh size, and that all the other eigen- 
values are bounded below and above nearly uniformly with respect to the jump and mesh 
size. As a result, we proved that the convergence rate of the preconditioned conjugate 
gradient method is uniform with respect to the large jump, and depends logarithmically 
on mesh size. These results ensure that the standard multilevel and domain decomposi- 
tion preconditioners are efficient and robust for finite element discretization of ( |1.1[ ) on 
quasi-uniform grids. In this paper, we extend our results to locally refined grids. 

The discontinuity of diffusion coefficients causes a lack of regularity of the solution 
to ( |1.1| ), which in turn, leads to deterioration in the rate of convergence for finite ele- 
ment approximations over quasi-uniform triangulations. Adaptive finite element meth- 
ods through local mesh refinement can be applied to recover the optimal rate of con- 
vergence [fTSll . In order to achieve optimal computational complexity in adaptive finite 
element methods, it is imperative to design fast algorithms for solving the linear system 
of equations arising from the finite element discretization. The distinct feature of apply- 
ing multigrid methods on locally refined meshes is that the number of nodes of nested 
meshes obtained by local refinements may not grow exponentially, violating one of the 
key properties of multilevel methods on uniform meshes that leads to optimal 0{N) 
complexity. Indeed, let N be the number of unknowns in the finest space, the complexity 
of multilevel methods with global smoothers can be as bad as 0{N'^) [|33ll . This prevents 
direct application of algorithms and theories developed in [i55il for quasi-uniform grids 
to locally refined grids. 

To achieve optimal 0{N) complexity, the smoothing step in each level must be re- 
stricted to the newly added unknowns and their neighbors (see [l6l[IIl|33l). Such methods 
are referred to as local multilevel methods in [6] . As an extreme case, one can preform the 
smoothing only on newly added nodes turning a coarse grid to a fine grid. The resulting 
method is known as the hierarchical basis method [|57l l8ll. In two dimensions, hierar- 
chical basis methods are proven to be robust for jump coefficient problems on locally 
refined meshes (cf. fSl). In three dimensions, however, classic multilevel and domain 
decomposition methods, including the hierarchical basis multigrid methods, deteriorate 
rapidly due to the presence of discontinuity of coefficients. To obtain robust rates of con- 
vergence for multigrid methods, one has to use special coarse spaces [|23l[39ll or assume 
that the distribution of diffusion coefficients satisfies the so called quasi-monotone con- 
dition [i23il . Therefore the three dimensional case is much more difficult. There are other 
works EII2EJ on optimal complexity of local multilevel methods in three dimensions, but 
the problems with discontinuous coefficients remain open. 

In this article, we shall design and prove the efficiency and robustness of local multi- 
level preconditioners for the finite element discretization of problem ( |1.1| ) on bisection 
grids - one class of locally refined grids. In these preconditioners, we use a global 
smoothing in the finest mesh; and for each newly added node, we perform smoothing 
only for three vertices - the new vertex and its two parents vertices (the vertices sharing 
the same edge with the new vertex). We analyze the eigenvalue distribution of the multi- 
level preconditioned matrix, and prove that there are only a fixed number of small eigen- 
values deteriorated by the coefficient and mesh-size; the other eigenvalues are bounded 
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nearly uniformly. Thus, the resulting preconditioned conjugate gradient algorithm con- 
verges uniformly with respect to the jump and logarithmically with respect to the mesh 
size of the discretization. We establish our results of this type in both two and three 
dimensions. 

To emploit the geometric structure of bisection grids, we use the decomposition of 
bisection grids developed in the recent work [[T9ll54ll . This approach enables us to intro- 
duce a natural decomposition of the finite element space into subspaces consisting only 
the newest vertices and their two parents vertices. In the analysis of these local multilevel 



preconditioners, one of the key ingredient is the stable decomposition (see Theorem 4.2 1. 
For the standard multilevel preconditioners on uniform mesh, in [|55ll we used the approx- 
imation and stability properties of the weighted projection (cf. W2\ ) to construct a 
stable decomposition. This weighted projection is no longer applicable for the local 
multilevel preconditioners, since it is a global projection. In order to preserve the local 
natural of the highly graded meshes, we introduce a local interpolation operator, which 



we manage to prove similar approximation and stability properties (see Theorem 3.4 and 
3.5) as the weighted L^-projection. Our local quasi-interpolation operator and the cor- 



responding analysis is more delicate than that in [fT9l |54] for the Poisson equation. We 
should remark that due to this space decomposition, we are able to remove the assump- 
tion, nested local refinement, which is used in most existing work on multilevel methods 
on local refinement grids [|2l |28l . 

The rest of the paper is organized as follows. In Section |2} we give some notation and 
recall some fundamental results as in ll55l . In Section]?} we study bisection grids, and 
review some technical tools from [fT9l[54l . Here we restrict ourself to a kind of special 
bisection scheme, namely the newest vertex bisection. Then in Section |4| we study some 
technical results of space decomposition, and present the optimal/stable decomposition 
and the strengthened Cauchy-Schwarz inequality on bisection grids. In Section[5} we an- 
alyze multilevel preconditioners, i.e., the BPX preconditioner and the multigrid V^-cycle 
preconditioner, and prove convergence results for the preconditioned conjugate gradient 
algorithm. In Section [6} we present numerical experiments to support our theoretical 
results. 

Throughout the article, we will use the following short notation, x <y means x < Cy, 
X > y means x > cy and x tz, y means cx < y < Cx where c and C are generic 
positive constants independent of the variables appearing in the inequalities and any other 
parameters related to mesh, space and coefficients. 

2. Preliminaries 

In this section, we introduce some notation, set up our problem, and review briefly 
some facts about the preconditioned conjugate gradient algorithm. 

2.1. Notation and Problem. Given a set of positive constants {ai}fi^, we define the 
following weighted inner products on the space H^{^1) 

M M 

{u, w)o,a = ^ ai{u, v)l-hq,), and (m, v)i^a = ^^(Vm, Vv)L2(n,) 

i=l i=l 

with the induced weighted norm || ■ ||o,a, and the weighted if^-seminorm | ■ 
respectively. We denote by 
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and the related inner product and the induced energy norm by 

(m, v)a = A{u, v) := (m, v)i^a, \\u\\a = \/A{u,u). 

To impose the Dirichlet boundary condition in ( |1.1[ ), we define 

^9dXd = {"^ ^ H^{il) : f Ir^ = fl'D in the trace sense}, 

and Hjj := Hq j.^ . Given a shape regular triangulation Th, which could be highly graded, 
we define Vh as the standard piecewise linear and global continuous finite element space 
on Th- Given / G H^^{Q) and G H^^'^{Tn), the linear finite element approximation 
of ( |1.1| ) is the function u G V/i fl H^^ , such that 

A{u,v) = {f,v)+ [ g^v, foraWv eVhnHh. (2.1) 

Given any mq G V/i fl iJ^^ p^, the problem ( |2.1[ ) is equivalent to finding m G V/i fl Hj^ 
such that 



v4(m,i;) = {f,v) + 



/ ^Tvt'- AK,t;), Vt; G V?,, n i/^. (2.2) 

We thus consider the space Vh,D '■= Vftflifl). The bilinear form A{-, ■) will then introduce 
a symmetric positive definite (with respect to standard L^-inner product) operator, still 
denoted by A, from Vh,D to Vh,D as 

{Au, v) = A{u, v). 

Define h G Vh,D as 

(6, t;) = (/, v) + I qnv - A{uq, v) Wv G Vh,D- 



N 



We then get the following operator equation on Vh,D 

Au = b. (2.3) 

For simplicity, in the remainder of the paper, we should omit the subscript D in Vh,D 
without ambiguity. 

We are interested in solving equation ( |2.3[ ) by the preconditioned conjugate gradient 
methods with BPX and multigrid preconditioners. Let us now review briefly some basic 
results concerning the preconditioned conjugate gradient method. 

2.2. Preconditioned Conjugate Gradient Method. Let i? be a symmetric positive def- 
inite (SPD) operator. Applying it to both sides of ( |2.3[ ), we get an equivalent equation 

BAu = Eh. (2.4) 



We apply the conjugate gradient method to solve ( |2.4[ ) and the resulting method is known 
as the preconditioned conjugate gradient (PCG) method, where B is called a precondi- 
tioned 

Let k{BA) = Amax(-BA)/Amin(-BA) bc the (generalized) condition number of the 
preconditioned system BA. Starting from an arbitrary initial guess uq, we have the fol- 
lowing well known convergence rate estimate for the fcth iteration Uk (k > 1) in PCG 
(see e.g. m) 



|m - MfclU ^ 2 ' V f^iBA) - 1 



I|m-moIU \^/k(ba) + i 

So if the condition number k{BA) is uniformly bounded, then PCG algorithm converges 
uniformly. Here the uniformity means the independence of the size of the matrix A. 
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Later on, when A is related to equation ( |1.1[ ), we shall also discuss the uniformity of 
convergence with respect to the jump of diffusion coefficients. 

If there are some isolated small or large eigenvalues, we can sharpen the above con- 
vergence rate estimate as stated in the following theorem. 

Theorem 2.1. [Sj Suppose that a{BA) = ao{BA) U ai{BA) such that there are m 
elements in ijQ{BA) and a < \ < (3 for each A G ai{BA). Then 



\u 



Uk A 



< 2K 



\u 



uo\\a 



where 



K 



max 

Ae(Ti{BA) 




k—m 



(2.5) 



/iecro(BA) 

If there are only m small eigenvalues in aQ^BA), say 



< Ai < As ■ ■ ■ < Ai < \m+i <■■■ <K 



then 



K 



n 

i=l 



A„ 



< 



An 



[k{BA) - 1)' 



(2.6) 



Therefore the convergence rate of PCG algorithm will be dominated by the factor ( ^J~j3ja- 
1)1 {^/JJa + 1), i.e. by (3 /a where (3 = A„(5A) and a = A„+i(5A). We define the 
"effective condition number" as follows. 

Definition 2.2. Let V be an n-dimensional Hilbert space and T : V — )■ V be a symmetric 
and positive definite operator. For any integer m E [1, n — 1], the mth effective condition 
number of T is defined by 

/rp\ Aniax(^) 

where \m+iiT) is the (m + l)-th minimal eigenvalue of T. 



As a corollary of Theorem 2.1 we have 

\\U - MfclU 



\u - Uo\\a 



< 2{k{BA) - ly 



'k„,{BA) - 1 

^Krn{BA) + 1 



k—m 



(2.7) 



From ( |2.7[ ), given a tolerance e, the number of iterations of the PCG method to reduce 
the relative error below the tolerance e is (cf. [|4l|5l) 



m + 



log 



+ m\\og{K{BA)-l)\ /co 



where cq = log (^{y^Km{BA) + 1) Km{BA) — 1)^ . Therefore if there exists an m > 

1 such that the mth effective condition number is bounded uniformly, then the PCG al- 
gorithm will still converge almost uniformly, even though the standard condition number 
K,{BA) might be large. 

To estimate the effective condition number, in particular \m+i{A), we use a funda- 
mental tool known as the Courant "minimax" principle (see e.g. Il24]| '). 
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Theorem 2.3. Let V be an n-dimensional Hilbert space with inner product (■, ■)v and 
T : V ^ V a symmetric positive operator on V. Suppose Ai < A2 < • ■ ■ < A„ are the 
eigenvalues ofT, then 

Xra+i[T) = max mm r — 

for z = 1, 2, ■ ■ • , n — 1. Especially, for any subspace Vq C V with dim(Vo) = n — m 

Am+i(T) > mm — . (2.8) 

If both A and B are SPD operators, then BA is SPD in the inner product induced by 
B~^ and A. Below, we shall apply Theorem 2.3 to T = BA and (m, f )v := {B~^u, v)l2. 
Therefore if we have an inequality of the type (Af , v) > c{B~^v, v) for all u in a suitable 
subspace Vq with dim(Vo) = n — m, we can get a lower bound of Xm+iiBA). 



3. Local Quasi-interpolation 

The theoretical justification of the robustness of multilevel preconditioners relies on 
establishing approximation and stability properties of certain interpolation operators. 
There are two difficulties: one is the locality and stability and another is the robustness 
with respect to the coefficient. 

The weighted L^-projection : L'^(Q) V/i defined by (Q^u, t;h)o,a = {u,Vh)o,a Vw^ e 
Vh was used in f 55ll59l for the case of uniform refinement. For the analysis of local mul- 
tilevel preconditioners, the interpolation operator should preserve certain local structure. 
Therefore, the weighted L^-projection, which is a global operator, is not appropriate. On 
the other hand, the standard nodal interpolation operator is local but not stable in the en- 
ergy norm. Local quasi-interpolation, such as Scott-Zhang operators [|411 . are developed 
to achieve both locality and stability. 

However, the stability constant will in general depend on the jump of diffusion co- 
efficients if we apply the standard quasi-interpolation globally on the whole domain. 
The value at a vertex is usually defined using a simplex in the patch of this vertex and 
thus depends on the diffusion coefficient in this simplex. For a vertex shared by several 
subdomains, this leads to the dependence of the ratio of coefficients. One remedy is to 
apply the quasi-interpolation on each subdomain and chose a sub- simplex in the quasi- 
interpolation. Indeed in the original paper [|4T1l . a {d — 1) sub-simplex is used. Such 
modification is suitable for the interior vertex relative to interfaces for which a common 
(d — 1) sub-simplex on the interface can be used to glue quasi-interpolations in different 
regions. For vertices on the boundary of the interface, i.e., edges in 3-D and vertices in 
2-D, in general there is no common {d — 1) sub-simplex but only {d — 2) sub-simplex. 
The trace of functions is not even well defined on {d — 2) sub-simplex. For exam- 
ple, the function value of a function at a point can be changed without changing this 
function. In the discrete level, it can be shown that the trace of a finite element function 
ona (d — 2) sub-simplex can be almost bounded by its Sobolev norm inside. Therefore 
we can simply set the function values at the vertices of {d — 2) sub-simplex to zero to 
glue quasi-interpolation operators defined in different domain. 

Below, we construct a quasi-interpolation operator by gluing Scott-Zhang operators in 
each subdomains and interfaces, and show that it is stable uniformly with respect to the 
jump of coefficients and nearly uniform to the mesh size of the triangulation. We stress 
that this local quasi-interpolation operator is designed for the analysis only, and is not 
needed in the practical implementation. 
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3.1. Notation on Triangulations. Let us introduce some notation related to the domain 
and its triangulations. As we mentioned earlier, we assume that the polygonal or poly- 
hedral subdomains Vli [i = 1, ■ ■ ■ ,M) are open, disjoint to each other, and satisfy 
^fii^i = ^- We denote Tij = dVli fl dVlj, or simply T if without ambiguity, as the 
interface between two subdomains fi, and Vtj. The subdomains Vti [i = 1, ■ • ■ M) may 
possibly have complicated geometry but we assume that they are resolved by an initial 
conforming triangulation %. Recall that a triangulation T is called conforming if the in- 
tersection of any two elements r and r' in T either consists of a common vertex, edge, 
face (when = 3), or empty. 

Let J\f, £ and J-' (when d = 3) denote the set of vertices, edges, and faces of T 
respectively. For each vertex p e Af,we define local patch Up := UrgpT and, for r E T, 
cUj- = UpgT-Wp. Similarly, on the (d — 1) dimensional interface F, Op, Og and oj denote the 
intersection of corresponding local patches and the interface. The linear finite element 
space associated to T is denoted by V(T), or simply V. More generally, for any subset 
S C T, V{S) denote the finite element subspace restricted to the subset G. Similarly, 
we should denote J^{G) C Af, S{G) C £ and -F(G') C as the set of vertices, edges, 
and faces in G C fi, respectively. 

For each element r E T, we define h^. = \t\^/^ and for the radius of its inscribed 
ball. In the whole paper, we assume that the triangulation is shape regular in the sense 
hr ~ Pt- Let h denote the piecewise constant mesh size function with h\r = hr, and 
/imin := minT-gT" hr- We should also denote he by the length of an edge e E £ and hfhy 
the diameter of a face f E T . Moreover, we define hp as the diameter of the local patch 
Up. By the shape regularity assumption, for all e, f,r C cUp, we have hp tz, h^ ^ hf tz, h.^. 

3.2. Technical Lemmas. For completeness here, we quote some technical lemmas from 
[fT2l . which will be used later for proving the approximation and stability of our local 
interpolation operator. 

In two dimensions, it is well known that H^{Vl) is not embedded into L°°(ri). But 
for finite element functions, we can control the norm by its H^-noxm with a factor 
I loe/) ■ 1 1/2 

I "mini 

Lemma 3.1 ( lfT2l Lemma 2.3]). For any subdomain Vti C M^, let V{VLi) be the finite 
element space based on a shape-regular triangulation T ofQi. Then for all v E V{Qi), 
it satisfies 



\v\\L°°{n,) < 



1 ^ 

log 



h 



mm 



1/2 



where Hi = diam(f2j) and h^m '■= min^-gT- hr- 

In three dimensions, the trace of an //^-function on an edge is not well defined. But 
for a finite element function, its L^-norm on an edge can be bounded by its if^-norm 



with a factor | log /imini It is a generalization of Lemma 3. 1 to three dimensions in the 
sense that controlhng the norm on a co-dimension 2 boundary manifolds. 

Lemma 3.2 ( [|T2l Lemma 2.4]). Given a polyhedral subdomain f2j C M^, let E cM. be 

any edge of Vti cind V(i7j) be a finite element space based on a shape-regular triangula- 
tion of Vti. Then for all v E V{Qi), there holds 

1/2 



Ml2{E) < 



1 

log 



^min 



where Hi = diam(fij 



8 



L. CHEN, M. HOLST, J. XU, AND Y. ZHU 



In the analysis of the local quasi-interpolation in Theorem |3 .41 below, we should ap 



ply Lemma 3.1 and Lemma 3.2 on each subdomain Qi, for which the diameter Hi 



diam(fij) ~ 1 is a fixed generic constant. 

3.3. Stable Local Quasi-Interpolation. Given a conforming triangulation % , the Scott- 
Zhang interpolation operator 11 : H^(^l) — )■ V{Th) can be defined as follows. For any 
p e J\f{Th), we choose a (rf — l)-simplex ap 3 p in Th- We remark that the choice of 
CTp is not unique (see Section [441 for the particular choice of for our purpose). Let 
{'^(7p,i : « = 1, ■ ■ ■ , (i} be the barycentric coordinates of Op. One can define the L^-dual 
basis {e„^^i : 2 = 1, ■ ■ ■ , c/} of {K^^i : i = 1, ■ ■ ■ , c?}, namely, /^^ Oaj„iKp,j = ^ij- We 
define a quasi-interpolation 11 as 



^^'= > . / 0,v\(Pp, (3.1) 




where {(l)p}peAf{Th) is the set of nodal basis of V{Th), and 9^^ = O^^^i. The following 
properties of the operator 11 can be found in [|4n [351 . 

Lemma 3.3. The interpolation operator 11 satisfies the following properties: 
(i) Stability: 

l|n^lli/i(r) < \\v\\h\u>^)] (3.2) 



(ii) Locality: 

(iii) Approximability: 



{Iiv)\r = v\r ifveViUr); (3.3) 



\h-\v-Uv)\\mr) < WvWmiu.^). (3.4) 



We apply the quasi-interpolation p.l[ ) on each subdomain, and denote Ilj : L^(fij) — )■ 



V(i7j) by the Scott-Zhang interpolation restricted to fi/. To be able to glue them together, 
we require for a vertex on the interior of the interface, we choose a common (c? — 1) sub- 
simplex shared by two sub-domains. By such choice. Hi and Uj will match on the vertex 
interior relative to the interface. 

We now define a local interpolation operator which has the desirable local approx- 
imation and stability properties in the weighted Sobolev norms. Given au E H^{^1), we 
define X> e V(7I) such that forp G 

rs,rr,V- / otherwise, 

For a vertex p, let ap be the (d — 1) -simplex chosen to define the nodal value at p. 
Then the interpolant is uniquely determined by the mapping p — )■ ap. In p.5[ ), if p is 
in the interior of some subdomain then ap C fij is chosen to be any {d — l)-simplex 
in T containing p; if p is in the interior of the interface T, then ap C T is chosen to be a 
{d— l)-simplex on the interface containing p. The choice of ap is not unique. However, in 
order to preserve the local structure of the adaptive grids, ap should be chosen carefully 
for each vertex p. This will be clear in Section |4] when we discuss the geometry of the 



bisection grids (see Section 4.4 for details). Now we are in the position to present the 
main result in this section: 
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Theorem 3.4. Let Vt dW^ with d = 2 or 3 and Tube a triangulation ofVt with mesh size 
h. Then for all u G H^{Q), we have 



Proof. Using the discrete Sobolev inequality Lemma 3J_ or 32_ on dT and the local H^- 
stability p.2| ) of Ilj, we have 



By the triangle inequality and the approximation property p.4| ) of Hi, we have 

\\h~\u-I^u)\\L2^n.^) 

< \\h-\u-U,u)\\mn,) + \\h-\UiU-I^u)\\mn,) 
^ \\u\\min,) + W^MlL^idT) 

TcdQi 

< \\u\\m{n,) + I log/lmin|5 • 

Multiplying by a suitable weight and summing up over all subdomains on both sides, 
we get the desired estimate. □ □ 

In general, we cannot replace „ by the energy norm \u\^^ in the above lemma; 
see [50 1 for a counter example. To be able to use ^ in the estimate, we introduce a 
subspace H}j{Q.) of H}j{Q) as follows: 

H}y{n) = lu e H}){n) : [ udx = for alH G / 1 , 

where / is the set of indices of all floating subdomains: 

I = {i: meas(afii n Yd) = 0}. 

Let mo := #/ be the cardinality of /. We emphasize that niQ is a constant, depending 
only on the distribution of the coefficients, and mo < M. In this subspace H}y(^l), the 
interpolation has the following properties. 

Theorem 3.5. For any v G H\){VL), we have the approximation property ofX'^ 

\\h-\v-Xlv)\\^^^<\logK,,,,[^\v\,^^, (3.6) 
and the stability ofX^ in the energy norm 

&ll,.< |log/lminN^;|i,,. (3.7) 

Proof. Voxv G Hjj^n), it satisfies the Poincare-Friedrichs inequality on each subdomain 



Therefore we get ||f ||o a ^ bli a • The inequality p.6[ ) then follows from Lemma 



3.4 



To prove inequality p.7[ ), we use the inequality p.6[ ) and the local projection Q^. : 



L^(r) — i- Vo{t) defined by Qru\r = \t\ ^ f^udx. Then on each element r G Th, we 
have 

< (ll^ -^>lli2(^) + \\V - QrV\\\2^^-^ ^ 
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where in the last inequahty, we used the approximation properties of Q^. Multiplying by 
a suitable weight and summing up over all r G T on both sides, we get 



where in the last step, we used inequality (3.6 1. □ □ 



Remark 3.6. When the coefficients satisfy the quasi- monotone assumption, the factor 
I log /iminl can be removed by arguments on a modified local patch; see [|23ll37l . □ 



4. Bisection Grids and Space Decomposition 

In this section, we give a short overview of the framework in the multilevel space 
decomposition on bisection grids in the recent work [[T9ll54ll . Most of the material in this 
section can be found there. 

4.1. Bisection Methods. We recall briefly the bisection algorithm for the mesh refine- 
ments. Detailed discussions can be found in ifTOl [1711331 and the references cited therein. 

Given a conforming triangulation T of ^l, for each element r E T,we. assign an edge 
of r to be the refinement edge of r, denoted by e(r) or simply e without ambiguity. 
This procedure is called labeling. Given a set of elements marked for refinement, the 
refinement procedure consists two steps: 

(1) bisect the marked element into two elements by connecting the middle point of 
the refinement edge to the vertices not contained in the refinement edge; 

(2) assign refinement edges for two new elements. 

Given a labeled initial grid To of ^2 and a bisection method, we define 

F(7o) = {T : T is refined from To by bisection method }, 
T(7^) = {T e F(7^) : r is conforming}. 

Namely F(7o) contains all triangulations obtained from To using the chosen bisection 
method. But a triangulation T G F(7o) could be non-conforming and thus we define 
T(7o) as a subset of F(7o) containing only conforming triangulations. 

Given any triangulation T, we define To = T, and the A;th uniform refinement 
Tk {k >1) being the triangulation obtained by bisecting all element in Tk-i only once. 
Note that for a conforming initial triangulation To with arbitrary labeling, Tk G IF(7o) 
but not necessarily in the set T(7o) in general. Throughout this paper, we shall consider 
bisection methods which satisfy the following two assumptions: 

(Bl) Shape Regularity: F(7o) is shape regular. 

(B2) Conformity of Uniform Refinement: Tk{%) G T{To) for all > 0. 

In two dimensions, newest vertex bisection with compatible initial labeling [|32l sat- 
isfies (Bl) and (B2). In three and higher dimensions, the bisection method by Kos- 
saczky [301 and Stevenson f43l will satisfy (Bl) and (B2). We note that to satisfy as- 
sumption (B2), the initial triangulation is modified by further refinement of each element, 
which deteriorates the shape regularity. Although (B2) imposes a severe restriction on 
the initial labeling, it is crucial to control the number of elements added in the comple- 
tion which is indispensable to establish the optimal complexity of adaptive finite element 
methods [[Ml. 
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4.2. Compatible Bisections. For a vertex p e M{T) or an edge e G £-{T), we define 
\hQ first ring of p or e to be 

7^p = {r G r|p G r}, 7^e = {r G T | e C r}, 

and the local patch of p or e as Up = UreTZp^, and = U^gT^^r. Note that tUp and are 
subsets of ri, while TZp and T^g are subsets of T which can be thought of as triangulations 
of Up and Wg, respectively. The cardinality of a set S" will be denoted by ^S. 

Given a labeled triangulation T, an edge e G i^(T) is called a compatible edge if e is 
the refinement edge of r for all r E TZe- For a compatible edge, the ring TZ^ is called a 
compatible ring, and the patch Wg is called a compatible patch. Let p be the midpoint of e 
and TZp be the ring of p in the refined triangulation. A compatible bisection is a mapping 
be ■ TZe ^ Tip- We then define the addition 

r + 6g ■.= r\neunp. 

For a compatible bisection sequence B := (&i, ■ ■ ■ , &fc), the addition T + i3 is defined as 

r+B = {{r+h) + b2) + --- + bk, 

whenever the addition is well defined. Note that if T is conforming, then T + &g is 
conforming for a compatible bisection be, whence compatible bisections preserve the 
conformity of triangulations. 

We now present a decomposition of meshes in T(7o) using compatible bisections. 



which will be instrumental later. We only give a pictorial demonstration in Fig. 4.1 to 
illustrate the decomposition. For the proof, we refer to ll54l . 



Theorem 4.1 (Decomposition of Bisection Grids). Let To be a conforming triangulation. 
Suppose the bisection method satisfies assumptions (B2), i.e., for all k > all uniform 
refinements Tk o/7o cire conforming. Then for any T G T(7o), there exists a compatible 
bisection sequence B = (&i, ■ ' ' j^n) with N = ^Af{T) — #A/'(7o) such that 

r = % + B. (4.1) 




Figure 4.1. A decomposition of a bisection grid. 

We point out that in practice it is not necessary to store B explicitly during the refine- 
ment procedure. Instead we can apply coarsening algorithms to find the decomposition. 
We refer to [20] (see also [fT8|| ) for a vertex-oriented coarsening algorithm and the appli- 
cation to multilevel preconditioners and multigrid methods. 

For a compatible bisection bi G B, we use the same subscript i to denote related 
quantities such as: 
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c * J • the patch of »j i.e. Wr, ; 
Cj: the refinement edge; « r ; . ^' V 

^, . , . ^ r- • Pi ,Pr '- two end points of 
Pi : the midpoint of Cj; ^t«'^n r 

0^ 0^. U Wp;^ U , ^ ^ diameter of w,-; 



7i = 7^) + ■ ■ ■ ,h 



1 5 



7?.j: the first ring of Pi in 7^. 



4.3. Generation of Compatible Bisections. The generation of each element in the ini- 
tial grid To is defined to be 0, and the generation of a child is 1 plus that of the father. The 
generation of an element r E T G IF (To) is denoted by g^. and coincides with the number 
of bisections needed to create r from Tq. For any vertex p E MiTo), the generation of 
p is defined as the minimal integer k such that p E MiTk) and is denoted by gp. In 
[|54l . we show that if 6j G i3 is a compatible bisection, then all elements of TZi have the 
same generation gi. Therefore we can introduce the concept of generation of compatible 
bisections. For a compatible bisection 6j : T^-g. — ?■ Tip., we define g^ = g^r), r E T^p.. 

Throughout this paper we always assume /;,(r) ~ 1 for r G To- Then since a bisection 
of a simplex will reduce the volume by half, we have the following important relation 
between generation and mesh size 

In particular, we introduce a "level" (or generation) constant L := laaXrergT- It is 
obvious that L ~ [| log /iminll ■ 

Different bisections with the same generation have disjoint local patches. Namely for 
two compatible bisections 6,; and bj with gj = gfj, we then have coi D tuj = 0. A simple 
but important consequence is that, for all u E L^(fi) and A; > 0, 



/i, ~7'% with 7= (^-j G(0,1). 



El 

9i=k 



^llo,a,aji ~ ll'"llo,a,o- (4-2) 



4.4. A Local Quasi-Interpolation. We define a sequence of quasi-interpolation oper- 
ators recursively. Let Xq : V(TAf) — ^ Vo be an arbitrary interpolation operator de- 
fined by p.5[ ). Assume Xf_]^ : V(TAr) — > V(Ti-i) is defined. Let bi be a compatible 
bisection, which introduces a new vertex pi from TI-i to % = T^-i + b^. We con- 
struct If : V(T/v) — > V(TI) as follows. If the new vertex pi E Td, we simply define 
{I°'v){pi) = to reflect the vanishing boundary condition of v. Otherwise, if pi ^ Yd 



we define the nodal value at pi through p.l[ ) with the choice of cTp. as follows: 

(i) if Pi is in the interior of some subdomain Qi, we choose a {d — l)-simplex (jp. 
containing p^; 

(ii) if Pi is in the interior of some interface F, we choose a (d — 1) -simplex ap. C F 
containing p^; 

(iii) otherwise, we simply let o"p. = and define (Xft>)(pj) = 0. 

For other vertices p G J\f{Ti^i), let (Xp G T^-i be the simplex used to define (Xf_^i;)(p), 
we update (X"t>)(p) according to the following two cases: 



(i) if cTp C uip{7i) we keep the nodal value, i.e., (Xff )(p) = (Xf_]^f )(p); 

(ii) otherwise we update dp as dp Wp(Ti) fl (Xp to define (Xff )(p). 



In either case, we ensure that the simplex cTp C ujp{Ti). In this way, we obtain a sequence 
of quasi-interpolation operators 

xf : v(r^)^v(T-), ^ = 0---Ar. 
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Note that in general X^t> 7^ v since the simplex used to define nodal values of Z'^v 
may not be in the finest mesh T/v but in T/v-i- Figure 14.21 illustrates the choice of (Xp in 
different cases in 2D. 




Figure 4.2. Update of nodal values Zfw to yield X°_^m: the element 
r chosen to perform the averaging that gives (XfM)(p) must belong to 
^pCTi)- This implies (Xf — Xf„^)M(p) 7^ possibly for p = Pi^Pi^^Pn and 
= otherwise. 



4.5. Stable Space Decomposition. Let G V(7i) denote the nodal basis at node 
p e Af{Ti). Motivated by the stable three-point wavelet construction by Stevenson H2l . 
we define the subspaces Vq = V(7o), and 

Let {(pp : p G A} be a basis of V{Tn), where A is the index set of the basis functions, 
and let Vp be the 1 -dimensional subspace spanned by the nodal bases associated to p in 
the finest grid. We choose the following space decomposition: 

N 

V:=^Vp + 5^V.. (4.3) 

Recall that bi only changes the local patches of two end points of the refinement edge Cj 
going from 7I-i to 71. By construction (Xf —I^_-^^)v{p) = forp G Af{Ti),p 7^ Pi,Pii or 
Pr^, which implies Vi := (X" —X^_^)v G Vj. Although I'^v 7^ t> in general, the difference 

V — X^t> is of high frequency in the finest mesh. Let us write v — X'^v = X^peA '-^^ 
basis decomposition. We then obtain a decomposition 

N 

V = ^Vp + ^Vi, Vie V^, (4.4) 

peA i=0 

where for convenience we define X":-^ := 0. Moreover, we introduce a subspace V := 

V n H}j{^l). Then we have the following stable decomposition. 

Theorem 4.2 (Stable Decomposition). Given a triangulation Tn = % + B in T(7o), let 
L = maXreTN fi'(^)- 



where Cd{L) 



14 L. CHEN, M. HOLST, J. XU, AND Y. ZHU 

(i) For any v E V, there exist Vp G Vp (p G A) and Vi E Vi {i = 1, ■ ■ ■ , N) such 
that V = Y.peA + Yj^=q 

N 

Yl K^lMla + HWla + Yl K'lHlla < Cd{L)\v\l,, (4.5) 
peA i=l 

L2, d = 2 
2^, d = 3 ' 

(ii) For any v E V, there exist Vp G Vp (p G A) and Vi E Vi {i = I, ■ ■ ■ , N) such 
that V = J2peA + T^iLo 

N 

Yl ^;'ll^pllo,a + ll^ollla + Y K'Mla < L^la (4-6) 
peA i=l 

Proof. The result of (i) is standard. We may use the standard nodal interpolation operator 
to define a decomposition using the hierarchical basis (cf. [[531 ). 

Now we prove (ii). Given a f G V, we define vq := X^v and Vi := {Xf — Xf_^)v. For 
v—X'^v = X^peA ^P' '^he approximability of the quasi-interpolation, cf. ( |3.6| ), we have 

E^'^il^^'llo.a ^ \\h-\v-X%v)\\l^ < L\v\l^. (4.7) 

pSA 



On the other hand, by Theorem 3.5 we obtain 



N 



i=l 

L 

< fE|l°g^-l) \\<a<L'\v\la- 



Then ( |4.6[ ) follows by adding the above inequality to inequality ( |4.7[ ). □ □ 



Remark 4.3. The estimate ( |4.5[ ) is not uniform for d > 2. For = 2, L ^ | log/imin| 
and the growth of C2{L) is acceptable. But for d = 3, the constant c^{L) = 2^ grows 
exponentially. This is the main reason that the hierarchical basis multilevel method de- 
teriorates rapidly in 3D (cf. [[581 [71[). For discontinuous coefficients problems, it seems 
unlikely to find a better decomposition with a better constants; see the counterexamples 

If the coefficients satisfy certain monotonicity, e.g. quasi-monotonicity (cf. [|23l [371) 
in the local patches, one can show that the interpolation operator defined above is stable 
in the energy norm without deterioration. □ 



Remark 4.4. With a close look at the proof of ( |4.6[ ), we may regroup the Vi = iXf—Xf_i)v 
into groups U^^G'(/) = {1, 2, ■ ■ ■ , N} such that for any i, j G G{l),ujj HUi = and 
therefore 

N L' 
E^'^^ll^^llo,'^-'^. = E E ^^^ll^j'llo.a,^. - ^'l^°S/imin||^^|?,a• 
The constant L' could be much smaller than L; see Section[6]for numerical examples. □ 
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4.6. Strengthened Cauchy-Schwarz inequality. An important tool in analysis of the 
multiplicative preconditioner is the following strengthened Cauchy-Schwarz inequality. 
A proof can be found in [fT9l[54ll . 

Lemma 4.5 (Strengthened Cauchy-Schwarz Inequality). For any Ui,Vi G Vj, i = 
0, 1, • • • ,N,we have 

N N 



J2 J2 

i=0 j=i+l 



N 



N 



< 



.1=0 



.1=0 



(4.8) 



As a corollary of ( |4.8[ ) and the inverse inequality, we have 



111, a 



h- \\ui 



■«IIO,a- 



i=0 



i=0 



(4.9) 



5. MULTILEVEL PRECONDITIONERS 

In this section, we shall analysis the eigenvalue distribution of the BPX preconditioner 
and the multigrid \^-cycle preconditioner on bisection grids, and prove the effective con- 
ditioner number is uniformly bounded. 

5.1. BPX (Additive) Preconditioner. To simplify the notation, we include Vat+i = V 
and rewrite our space decomposition as V = J2fJo^ ^i- Based on this space decomposi- 
tion, we choose SPD smoothers Ri : Vi — )■ V, satisfying 

{R^^Ui,Ui)o^a ^ hr^{ui,Ui)o,a, Vwj G Vj. (5.1) 

According to [55], both of the standard Jacobi and symmetric Gauss-Seidel smoother 
satisfy the above assumption. On the coarsest level, i.e. when i = 0, we choose the exact 
solver i?o = ^. Let : V — Vi be the weighted projection. Then we can define 
the BPX-type preconditioner 

N+l 
i=0 

It is well known [|49l EH HH that the operator B defined by (|5]2]) is SPD, and 

Af+l 

{B~'^v,v)o^a = inf y^(R~'^Vi,Vi)o^a- (5.3) 

We have the following main result for the BPX preconditioner. 

Theorem 5.1. Given a triangulation Tn = To + B in T(7o), let L = max-r^TN di'^)- ^^^^ 
the BPX preconditioner defined in ( |5.2[ ), we have 

k{BA) < CiCd{L), andnmo{BA) < CqL^ . 

Consequently, we have the following convergence estimation of the BPX preconditioned 
conjugate gradient method: 

II II //^ r 1 \ k~mo 



\u-uo\\a \CoL 
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Proof. First of all, let us estimate Amax(-B^)- For any decomposition v = v+^-^q Vi, v G 
V,Vi G Vi, we have 

AT „ N N+1 



|2 < 



I„~,ll2 



i=0 



i=0 



i=0 



In the second step, we used the inverse inequality and the inequality ( |4.9[ ). In the third 



step, we used the assumption ( |5.1[ ) of Ri . Taking infimum, we get 

Af+l 



V,V On, 



which implies that 



EiV 

■XBA) < 1. 



i=0 



To estimate Amin, in view of (5.3) we choose the decomposition as in the stable de- 



composition Theorem 4.2 (see ( |4.5[ )) to conclude that 

7V+1 

{B'^V,v)o^a < '^{Ri^Vi,Vi)o^a < Cd{L){Av , v)o,a, 
i=0 

which implies that X^in{BA) > Cd{L). Therefore we have k{BA) < Cd{L). 

On the other hand, if we app ly (|4.6[ ) in the subspace V C V, we obtain Xm^^^i{BA) 



> 



L by the "min-max" Theorem 2.3 Hence we get an estimate of the effective condition 



number Kmo(-BA) < L^. The convergence rate estimate then follows by Theorem 2.1 
This completes the proof. □ □ 

From this convergence result, we can see that the convergence rate will deteriorate 
a little bit by q(L) as L grows. But since rriQ is a fixed number, when k grows, the 
convergence rate will be controlled by the effective condition number, which is bounded 
uniformly with respect to the coefficient and logarithmically with respect to the mesh 
size. Notice that L ~ | log /imin| and thus the asymptotic convergence rate of the PCG 



algorithm is 1 — 



1 



C| log hniinl 



for h <l. 



Remark 5.2. The estimate k{BA) < CiCd{L) is sharp in the sense that there exists an 
example on BPX preconditioner such that k{BA) ~ (cf. [i36J). □ 

Remark 5.3. Here we should emphasize that the convergence rate estimate in Theo- 
rem |5.1| holds for general substructures. In some special circumstance, for example 



"edge type" or "exceptional" in the terminology in [|36ll . or "quasi-monotone" coeffi- 
cient in [23 J, we can sharpen the convergence estimate in Theorem |5.1| by a modification 
of Theorem |4.2[ see [36|. □ 

5.2. Multigrid (Multiplicative) Preconditioner. We shall use the following symmetric 
V-cycle multigrid as a preconditioner in the PCG method and prove the efficiency of such 
a method. Let Ai := A\\;-. Then one step of the standard V-cycle multigrid B : V V 
is recursively defined as follows: 

Let Bq = Aq^, for i > and g G Vj, define B^g = w^. 

(i) Presmoothing : wi = Rig; 

(ii) Correction: ^2 = ^1 + -Bi-iQi-i(fi' - AiWi); 

(iii) Postsmoothing: ^3 = ^2 + R*{g — AiW2). 
Set B = Bn+i- 



For simplicity, we focus on the case of exact subspace solver, i.e., Ri = A^^ for 
i = 0, ■ ■ ■ , and for the finest level, Rn+i is chosen as Gauss-Seidel smoother, which 
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can be also understood as the multiplicative method with exact local solvers applied to 
the nodal decomposition [|5T1l . Let Pp : V — )■ Vp and Pj : V — )• be the orthogonal 
projection with respect to the inner product (■, ■)«. For our special choices of smoothers, 
we then have 



p6A 



N 



N 



BnA 



\I -BA\ 



,1=0 

N 



vi=0 



pGA 



i=0 



For exact local solvers, we can apply the crucial X-Z identity [|56l to conclude 

1 



\I -BA\ 



1 - 



1 + co' 



(5.4) 



where 



Co = sup inf 

II''IIa=1 ^=EpeA -"P+E^o ' 



N N 

j=i+i p£A peA 



J=0 



q>p 



Theorem 5.4. Given a triangulation Tn = % + B in T(7o), let L = max^-gT-^ 9{'^)- For 
the multigrid V -cycle preconditioner B, we have 

k{BA)<c,{L), n^,{BA)<L\ 

Consequently, we have the following the convergence rate estimate of the BPX precondi- 
tioned conjugate gradient method: 



\U - Uk\\A 
\U - MolU 



< 2(Ciq(L)-1) 



mo 



CqL-I 
CoL + 1 



Proof. Since / — BA is a non-expansive operator, we conclude Amax(-BA) < 1. Since 
/ — BA is SPD in the A-inner product and Ainax(-BA) < 1, we have 

\\I-BA\\a = max{|l-A^in(PA)|,|l-A^ax(5A)|} = l-A^in(Pv4). 

To get an estimate on the minimum eigenvalue of BA, we only need to get a upper bound 
of the constant cq in (|5.4[). 



To do so, for any G V, we chose the decomposition in Theorem 4.2 That is. 



N 



V = v + ^Vi, with vo = I^v, Vi = (Xf - Xf_ Ju, 

1=1 

where v = v — X%v = J2peA ^p- Then by shape regularity of the triangulation, we have 



Co < 



N N ^ N 

i=0 i=i+l i=0 peA 



q>p 
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We estimate these three terms as follows. For the last term, by the finite overlapping of 
nodal bases, we have 



i\\a ~ 



EIIE 



l\\A,t 



peA 



q>p 



peA q>p 

< \^ \\v l|2 < /7-2|L, ||2 

~ /_^\\^p\\A,uip ~ / J V II "^P Ho, a,tjp 



peA 



peA 



< \\hr\v -xiv)\l^<\\v\ 



For the middle term, we regroup by generations and use (4.2) to get 



N 



E k-' 



E E i^'* 

fc=0 l,gi=k 
L 



< 



E 

k=0 



< 



„~,l|2 



EE 

k=0 l,gi=k 



Ml., 



Ma = L\Ma- 



For the first term, we define = Pi ^j) ^^'^ := -Po(^ — ^o) and apply the 

strengthened Cauchy Schwarz inequality, cf. Lemma [43] to get 



N 



N 



Ek^E 



i=0 



j=i+i 



|2 

\a 



N N 
i=0 j=i+l 

N 



< 



-2|| ||2 



i=l 

< c,{L)\\v\\l 

Here the constant Cd{L) can be improved to if we consider the decomposition ( |4.6| ) of 



V eV. Combined with the Mini-Max Theorem 



2.3 



yields 



and thus 



k{BA)<c,{L), k^,{BA)<L\ 



Finally, the convergence rate of the PCG method follows by Theorem 2.1[ □ 



□ 



Follow the same proof as Theorem 5.4 we can also obtain the following convergence 
result for the local multigrid y-cycle solver. 

Corollary 5.5. For the multigrid V -cycle algorithm defined above on bisection grids, we 
have 

\\E\\a = \\I-BA\\a = 1-^^. 



Co 



where Cq < 



This corollary implies that multigrid alone is not robust, especially in 3D. In this case, 
the convergence rate of multigrid will be proportional to 1 — 2^^ c:^ \ — h^^^, which 
deteriorates rapidly as the mesh size become small. Remark 5.3 is also applicable here, 
i.e., all the above estimates are estimates for the worst case. For the special circumstances 
mentioned in Remark [53| the estimates can be improved in the same way. 
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6. Numerical Experiments 

In this section, we present some numerical experiments to support the theoretical re- 
sults in previous sections. In the implementation of the adaptive loop, we use a modifica- 
tion of the error indicator presented in [1371 . Some other a posteriori error indicators for 
jump coefficients problem ( |1.1[ ) can be found in ||9l|2Tl|44l[i4i|. The adaptive algorithm 
using different error indicators will generate different grids. However, we emphasize that 
the robustness of the local adaptive multilevel preconditioners is independent of how the 
grids are generated in the refinement procedure. 

The implementation of the BPX preconditioner and the multigrid methods are stan- 
dard, and can be found in, for example, [fT3l |53,l . The implementation of the PCG al- 
gorithm can be found in [|24l l38l . All numerical examples are implemented by using 
iFEM [fTSl . We only present three-dimensional examples here and refer to [|20l for two- 
dimensional ones. In the PCG algorithm, we use the stopping criterion 

II h i"— 1 II 



In the implementation of the local multilevel preconditioners, we use an algorithm for 
coarsening bisection grids introduced by [20] for two dimensional case and [ 1 8] for three 
dimensional one. The coarsening algorithm will find all compatible bisections and re- 
group them, with possibly different generations, into groups U[^^G{1) = {1, 2, ■ ■ ■ , N} 
such that for any i,j E G{l),ujj D Ui = 0. Each coarsening step is corresponding to a 
level in the multilevel terminology, and the total number of levels is L'. There are two 
major benefits of using this coarsening algorithm. 

(i) We do not need to store the complex bisection tree structure of the refinement 
procedure explicitly in the algorithm. Instead, we only need the grid information 
on the finest level and the coarsening subroutine will restore multilevel structure. 

(ii) Our numerical evidence shows that the number of nodes will decrease around 
one half in one coarsening step. Therefore the constant L' is much smaller than 
the maximal generation L ~ | log \ ■ 

In what follows, we will use some shorthand notation for the different algorithms imple- 
mented. 

• TPSMG stands for the V-cycle multigrid with Three-Point Smoothing (TPS), 
which only performs smoothing on new vertices and their two direct neighbors 
sharing the same edge. 

• TPSMGCG is the PCG algorithm using the TPSMG as preconditioner. 

• TPSBPXCG is the additive version of TPSMG preconditioner. 

Among all these algorithms, the main focus of this paper is the behavior of TPSMGCG 
and TPSBPXCG. In the numerical experiments below, we also report some results for 
TPSMG for comparison. 

Inspired by [|36l [50l [55l . we consider solving the model equation ( |1.1[ ) in the cubic 
domain f2 = (—1, 1)^. Let the coefficient a(x) be the constants ai = 02 = 1 and = e 



on the three regions Qi, Q2 and Q3 respectively (see Figure [6T| ), where 

n-^ = (-0.5, 0)^ n2 = (0, 0.5)^ and ^3 = \ (Hi U Ha). 
We choose / = 1 and impose the following boundary conditions: Dirichlet conditions 

M{-l}x[-l,l]x[-l,l] = O5 ""{Ijxl-l.llxl-l,!] = 

and homogenous Neumann boundary conditions on the remaining boundary. For this 



problem, singularities occur along edges of fii and ^72• Figure 6.2 shows an adaptive 
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Figure 6.1. The coefficients ai = a2 = 1 in the gray domains Vti and 
^2, and as = £ in the rest of the domain. 

mesh and the corresponding finite element approximation after several iterations of the 
adaptive algorithm. To view the mesh around the singularity, we only show half of the 
domain VL. 




I 



4500 
4000 



Figure 6.2. An adaptive mesh and finite element solution with e 
and 36466 vertices. 



Tables 6.1 ■ 6.4 give comparisons of the number of iterations for three different algo- 



rithms: TPSMG, TPSMGCG and TPSBPXCG algorithms, respectively, with the choice 
of e = 10""^, 10~^, 10^ and 10^. As we observe from these tables, the number of itera- 
tions for TPSMG algorithm grows rapidly as the mesh is refined when e is small. On the 
other hand, the number of iterations for TPSMGCG and TPSBPXCG is very robust and 
only grows a little bit when the mesh is refined, as we expected from the theory. We also 
observe that if e is large, the TPSMG algorithm will converge uniformly. This is because 
the coefficient in Vt^, which contains the Dirichlet boundary, is dominant. In this case, 
we could use the standard multigrid analysis (as in [51]) to show the robustness of the 
preconditioners. 



Figure [63] shows the eigenvalue distributions for the TPSMGCG and TPSBPXCG 
preconditioned systems. As we can see from the figure, there is one small eigenvalue 
for both preconditioned systems. This agrees with the theoretical results, the number of 
small eigenvalues is bounded by the number of floating subdomains niQ = 2. 



Figure 6.4 shows the condition number and effective condition number of TPSBPXCG 
and TPSMGCG preconditioned systems. From Figure |6.4[ we observed that when e 
is small, the condition number deteriorates (k(BA) E [3, 1100] for TPSBPXCG, and 
k{BA) e [3, 125] for TPSMGCG as we can see from the figure). On the other hand, if 
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DOF 


TPSMG 


TPSMGCG 


TPSBPXCG 


4913 


41 


13 


18 


5505 


62 


15 


18 


6617 


89 


18 


21 


8666 


99 


19 


19 


10585 


98 


19 


20 


12411 


125 


23 


25 


16353 


154 


23 


23 


21248 


182 


22 


23 


27755 


197 


26 


32 


36466 


178 


27 


29 


43271 


238 


25 


30 


51163 


283 


28 


36 


72349 


395 


32 


34 


89146 


424 


31 


34 


104747 


413 


34 


38 



Table 6.1. Comparison of Number of Iterations for TPSMG, 
TPSMGCG and TPSBPXCG when e = 10 ^ 



DOF 


TPSMG 


TPSMGCG 


TPSBPXCG 


4913 


46 


13 


17 


5550 


51 


15 


17 


6743 


61 


17 


20 


8907 


65 


16 


19 


10729 


66 


17 


20 


13281 


86 


20 


24 


17146 


90 


20 


21 


23139 


90 


20 


24 


28613 


160 


25 


29 


37338 


175 


24 


27 


43610 


149 


22 


26 


52715 


154 


25 


31 


72967 


238 


28 


29 


89320 


165 


25 


33 


113131 


294 


30 


38 



Table 6.2. Comparison of Number of Iterations for TPSMG, 
TPSMGCG and TPSBPXCG when e = IQ-^. 



we get rid of the first small eigenvalue, the effective condition number ki{BA) of TPS- 
BPXCG and TPSMGCG preconditioned systems (the black and red lines, respectively) 
are almost identical for different e. This indicates that the effective condition numbers 
are uniform with respect to the jumps. Moreover, as we can see from Figure 6.4 ki (BA) 
are mildly increasing with respect to the DOFs G [1, 80] for TPSBPXCG, and 

Ki{BA) E [1, 30] for TPSMGCG). These results agree with our theoretical expectations 
from Section |5l 
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DOF 


TPSMG 


TPSMGCG 


TPSBPXCG 


4913 


16 


10 


14 


5279 


37 


15 


15 


5867 


43 


17 


18 


6522 


48 


16 


19 


7562 


68 


17 


18 


9493 


61 


17 


18 


11858 


49 


15 


18 


15257 


68 


15 


18 


20649 


61 


16 


19 


27946 


49 


17 


21 


36735 


52 


16 


20 


48890 


58 


16 


22 


68297 


71 


18 


22 


89872 


55 


16 


21 


119109 


61 


17 


23 



Table 6.3. Comparison of Number of Iterations for TPSMG, 
TPSMGCG and TPSBPXCG when e = 10^. 



DOF 


TPSMG 


TPSMGCG 


TPSBPXCG 


4913 


16 


10 


14 


5269 


37 


15 


15 


5863 


42 


17 


18 


6493 


45 


16 


18 


7531 


68 


17 


18 


9419 


59 


16 


17 


11721 


46 


15 


18 


14941 


69 


15 


18 


20065 


59 


16 


19 


27199 


47 


17 


21 


35601 


59 


16 


20 


47743 


55 


16 


22 


66989 


71 


18 


21 


88079 


57 


16 


21 


116739 


56 


17 


23 



Table 6.4. Comparison of Number of Iterations for TPSMG, 
TPSMGCG and TPSBPXCG when e = 10^ 



7. Conclusion 

In this paper, we designed local multilevel preconditioners based on the decomposition 
of the finite element space into 3-point subspaces for the highly graded mesh obtained 
from adaptive bisection algorithms. To analyze the behavior of the local multilevel pre- 
conditioners, we introduced a local interpolation operator and proved some approxima- 
tion and stability properties of it. Based on these properties, we showed the decompo- 
sition of the finite element space is stable, which is a key ingredient in the multilevel 
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10 15 20 25 



(a) Eigenvalues for TPSBPXCG 



(b) Eigenvalues for TPSMGCG 



Figure 6.3. Example 2: Eigenvalues of BA when e = 10 ^ with 12411 vertices 



(Effective) Condition Number for BPX preconditioner 



(Effective) Condition Number for IVIG preconditioner 



-k(BA), s =10"' 
k(BA), s =10"^ 
_K,(BA), s=10" 

_K^(BA), s=10" 



i^****" 





Degree of freedoms 



(a) TPSBPXCG 



2.5 
X 10' 




Degree of freedoms 



(b) TPSMGCG 



Figure 6.4. Example 2: k,{BA) and k,i{BA) for the cases e 
10-^ 10"^ w.r.t the DOFs. 



analysis. This enabled us to analyze the eigenvalue distributions of the preconditioned 
systems. In particular, we showed that there are only a small fixed number of eigenval- 
ues that are deteriorated by the coefficients and mesh size, and the other eigenvalues are 
uniformly bounded with respect to the coefficients and logarithmically depends on the 
mesh size. As a result, we proved the asymptotic convergence rate of the PCG algorithm 
is uniform with respect to the coefficient and nearly uniform with respect to the mesh 
size. Moreover, the overall computation complexity of these multilevel preconditioner 
are nearly optimal. Numerical experiments justified our theoretical results. 
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